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We study the stability of vortex-lines in trapped dilute 
gases subject to rotation. We solve numerically both the 
Gross-Pitaevskii and the Bogoliubov equations for a 3d con- 
densate in spherically and cilyndrically symmetric stationary 
traps, from small to very large nonlinearities. In the sta- 
tionary case it is found that the vortex states with unit and 
m — 2 charge are energetically unstable. In the rotating trap 
it is found that this energetic instability may only be sup- 
pressed for the m = 1 vortex-line, and that the multicharged 
vortices are never a local minimum of the energy functional, 
which implies that the absolute minimum of the energy is not 
an eigenstate of the L z operator, when the angular speed is 
above a certain value, > Q2 ■ 

PACS: 03.75.-b, 02.70.Hm, 03.65.Ge 



I. INTRODUCTION 

Since the first experimental realization of Bose- 
Einstein condensation (BEC) in weakly interacting gases 
W , there has been a huge theoretical and experimental ef- 
fort to study its properties in the framework of fully quan- 
tum theories and in the so called mean field limit (Gross- 
Pitaevskii -GP- equations) . These equations are formally 
Nonlinear Schrodinger Equations (NLS) ||] which appear 
in many fields of physics, e.g. in bulk superfluids and 
nonlinear optics to cite only a few examples. 

All of these physical systems have been long known 
to exhibit solutions corresponding to topological defects 
P,PJ, one of the simplest being known as vortices (in 
two spatial dimensions) or vortex-lines (in three spatial 
dimensions). Vortices are localized phase singularities 
with integer topological charge which analogous in the 
hydrodynamic interpretation to vortices those appearing 
in fluid dynamics j|] . In the framework of BEC studies it 
has been raised the question of whether these nonuniform 
clouds of condensed gases may support the existence of 
vortices in a stable form. 

There is a huge literature on vortices and vortex prop- 
erties in the framework of NLS equations (including its 
particular cubic version, the GP equation) and its non- 
conservative extensions, the Ginzburg-Landau (GL) sys- 
tem, and vector GL models. In particular the stability 
of m-charged GP vortices in two dimensions was studied 
in In three dimensions the GL case has been re- 
cently considered and vortex lines geometric instabili- 
ties have been found to strongly deform the vortex lines. 
However the GP equation cannot be obtained as a limit 



of the GL studied there since dissipation and diffusion 
are essential ingredients of the models studied in Ref. . 
This fact makes the conservative case (GP) interesting 
by itself. Other analysis of vortices and vortex stability 
in the framework of Nonlinear Optics are included in Ref. 

The current setups utilized to generate Bose-Einstein 
condensates use a magnetic trap to confine the atomic 
cloud which is modeled by a parabolic trapping potential. 
This is a distinctive feature from the common NLS sys- 
tems, in which the vortices are free and move in an homo- 
geneous background. The dynamics of a vortex in a spa- 
tially inhomogeneous two dimensional GP problem was 
studied in Ref. § using the method of matched asymp- 
totic expansions, but the authors did not consider the 
stability of the 2D vortex itself. In principle, the vortex 
motion equations [|| can be used to study the motion 
of a single 2D point vortex in spatially inhomogeneous 
GP problems. However, the dynamics of the many vor- 
tex case is more complicated and by no means trivial. 
For simple approaches to the problem which do not in- 
clude the effect of vortex cores on the background field 
see Ref. |k]]. The dynamics of 3D vortices is yet more 
complicated allowing the so-called reconnection. To our 
knowledge there are no analytical results but only quali- 
tative numerical observations available |ll[] . Another the- 
oretical framework where non-homogeneous dynamics of 
vortices has been investigated is the possibility of pin- 
ning vortices in type-II superconductors ]l2| ], but it is 
only the dynamics has been considered through analyt- 
ical approximation techniques with no comparison with 
numerics. In all the previously discussed cases the vortex 
stability is given for granted. 

In the framework of Bose-Einstein condensed gases 
studies the problem of the vortex stability has been con- 
sidered in various papers that try to solve the problem of 
linear and global stability, either from a purely analytical 
point of view, such as in |13|-|l5| , or by mixing analytic 
and numerical techniques, |17| , |lq |. 

In Ref. |l7j the authors solve the Gross-Pitaevskii 
equation and find the energies of the condensate in vortex 
states, for a number of particles up to N = 10 4 . In Ref. 
Jl8| the authors solve the Bogoliubov equations for an 
unit charge vortex in a stationary trap with axial sym- 
metry, their results being also limited to N < 10 4 . In 
Ref. jl9| the authors address the problem of minimizing 
the energy functional with a reduced basis of trial states 
that is only valid in the limit of small U. 

In this paper we unify and substantially extend what 
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has been done in previous works regarding these two 
questions: global energetic stability and local stability of 
vortex states. First, in Sect. |l| we solve the GPE for an 
axially symmetric harmonic potential, with or without 
the action of a uniform magnetic field which resembles 
the effect of a rotating trap. We calculate the lowest sta- 
tionary solutions that have a well defined value of the 
third component of the angular momentum, m = (L z ), 
and we do this for small and for very large values of the 
nonlinearity (N ~ 10 7 ). We find that there are contin- 
uous intervals of the "angular velocity", |T2 m ,f2 m+1 ), in 
which the m-charged vortex state becomes energetically 
stable with respect to other states of well defined vortic- 
ity. In Sect. [II, we study Bogoliubov's equations from 
two different points of view: as a consequence of a linear 
stability analysis of the Gross-Pitaevskii equation (GPE), 
and as the first corrections to the mean field theory of the 
dilute condensate. The concepts of dynamical and ener- 
getic stability are defined, and it is demonstrated that 
any possible destabilization of the system must be either 
of energetic nature, or grow polynomially with respect to 
time. We next solve the Bogoliubov equations for m = 1 
and m = 2 unperturbed vortex states in stationary traps. 
It is found that the m = 1 and m = 2 vortices are only 
energetically unstable, which means that the lifetime of 
both configurations is only limited by dissipation. A sim- 
ilar treatment reveals that rotation can only stabilize the 
unit charge vortex-line if the angular speed is in a suit- 
able range, £1 £ [Six, ^2)5 while outside of this range, 
il 2 < < fl c , the minimum of the energy functional is 
not an eigenstate of the L z operator -i.e, it is not sym- 
metric under rotations. These results are confirmed by 
numerical simulations of the evolution of perturbed vor- 
tices. In Sect. |^ we summarize our work and discuss 
their implications. 



\tl>\ J d?x = 1. 



(2) 



It is convenient to express Eq. ([!]) in a natural set of 
units which for our problem is built up from two scales: 
the trap size (measured by the width of the linear ground 
state), ao = \/h/mu>, and its period, r = l/u>. With 
these definitions the equation simplifies to 
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while maintaining the normalization. 

The new parameters, fl = Ml and U = AnNa/ao, rep- 
resent the "angular speed" of the trap and the adimen- 
sionalized interaction strength, respectively. For stability 
reasons (see below) , f2 will be of the order of magnitude 
of or smaller than the strength of the trapping, u. The 
other parameter, U, will take values from to 6 x 10 4 . 
As of the experiments with rubidium and sodium, this 
implies a minimum of 10 6 and a maximum of 10 7 atoms 
which is in the range of current and projected experi- 
ments. The shape of the trap is dictated by the geometry 
factor, and in this work it will typically take two possible 
values: 7 = 1, corresponding to a spherically symmetric 
trap, and 7 = 2, corresponding to an axially symmetric, 
elongated trap. 

A stationary solution of (|J) will be of the form 



ip(x, t) 



(x) |2C|] , where fi may be interpreted 



both as a frequency and the chemical potential 



iA + ^ + i(7V + z 2 ) + [/|0| 2 



4>. (4) 



Any solution of (Q) has an energy per particle which is 
given by the functional 



II. VORTEX SOLUTIONS OF THE GPE 

A. Stationary states of GPE in a uniform and 
constant magnetic field 

For small temperatures and small densities, the con- 
densate is modeled by the Gross-Pitaevskii equation 
(GPE) We will always refer to an axially symmet- 

ric trap with term that accounts for rotation around the 
Z axis and may be generated by a weak magnetic field. 
The form of the equation is 



ih 
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dt 2m T ' 2 
+ U N\if)\ 2 il> + QL z i) 
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(i) 



Here Uq = Airh^a/m characterizes the interaction and 
is defined in terms of the ground state scattering length 
a. In all cases we will take the normalization condition 
to be 



E(r.A') = / ( i|VVf-^cW 



7 v+z 2 + c/|v>r) iv-r 



For a stationary solution it becomes 

£(V,iV)= M -f / I0I 4 



(5) 



(6) 



The stationary solutions of (||) may also be interpreted 
as the minimization of 



(7) 



subject to the constrain of Eq. (j2j). In that case /j, is 
nothing else but the Lagrange multiplier of the norm. 

Since we are interested in single vortex solutions to 
the GPE we will restrict our analysis to stationary states 
that are also eigenstates of the L z operator. That is, 
we will look for solutions of the form ip(r,z,9,t) — 
e~ lfJ - t e lme <f>(r, z). Summarizing, our goal will be to find 
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the unit norm functions (j) ™' (r, z) and real numbers fi 
which are solutions of the equation 



iA-mQ+i( 7 2 r 2 + * 2 ) + ^)| 2 



(8) 



Our treatment on the following sections will be fully 
three-dimensional and no spurious conditions (e.g. pe- 
riodicity) will be imposed on the boundaries. We want 
to obtain at least the lowest energy state for each value of 
the vorticity, m. Also the dependency of spectrum with 
the nonlinearity and the angular velocity, f2, is interest- 
ing since it will allow us to find whether the vortex-line 
states may become energetically favorable. 



B. Numerical method 

Due to the nonlinear nature of the problem we want to 
solve [Eq. (||)] there are not many analytical tools avail- 
able. The most common (and maybe easiest) approach 
to the problem is to discretize the spatial part and per- 
form time evolution in imaginary time while trying to 
preserve the normalization, a method which is related to 
the steepest descent. The precision of the solution de- 
pends of the type of spatial discretization used: finite 
differences (used for example in Refs. [[HJpT]]) or spectral 
methods (such as the one used in Ref. p3|). However 
these common methods, such as finite differences [fl7| and 
similar spectral methods |l8| ], have reached a maximum 
value of the interaction of U — 10 3 , which should be 
contrasted with the value U = 10 5 that can be attained 
which the technique to be presented later. 

Properly speaking our technique is a Galcrkin type 
method where one performs the expansion of the un- 
known solution in a complete basis of the Hubert space 
under consideration. For convenience we have used the 
basis of eigenstates of the harmonic oscillator with fixed 
vorticity. In this basis our stationary solution is ex- 
pressed as 

tfW(af, t) = e-^e ime £ c n P^\r, z), (9) 

n 

Here the single index, n, denotes two quantum numbers, 
(n z ,n r ), regarding the axial and radial degrees of free- 
dom, and P« is a product of a Hermite polynomial, a 
Laguerre polynomial and a Gaussian 



P^ = C n H nz {z)L nr {p 2 )r m e ime e-^ 2+ ^'\ (11) 
C = 
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v /7y / 7r2" z n z ! y tt (n r + m)! 



(12) 



with p = rj yfy 

Next, following the same convention about the indices, 
we have introduced this expansion into Eq. (ph to obtain 



(e*° -nm-^ct + uY, 42te<*<* = o» ( 13 ) 
jkl 

where E\° is the harmonic oscillator energy of the mode 
Pi m and the tensor A l ^ kl has the following definition 



(m) p(m) p (m) p (m) 



(14) 



bmce the P i (m) are products of known polynomials by 
exponentials it could be possible, in principle, to evalu- 
ate the coefficients in the tensor exactly with a Gaussian 
quadrature formula of the appropriate order. This ap- 
proach was used in Ref. |16[ for the three-dimensional 
case. However, when one wishes to use a large number of 
modes (which in our case is of about 1600 for each value 
of m) to achieve large nonlinearities, the search of the 
quadrature points becomes more difficult than perform- 
ing a stable integration by means of some other methods, 
of which the simplest accurate one is Simpson's rule [g2| . 

Once we fix all of the constants, E ho , A\j^, [i and 
a guess for the solution, it is feasible to solve ( p"3| ) it- 
eratively -e.g. by Newton's method p2|j . However, it is 
wiser to perform two simplifications before implementing 
the algorithm. The first one is that all of the eigenfunc- 
tions, Pn > can be made real and thus we can impose 
the coefficients in the expansion, {c„}, to also be real. 

The second optimization is that, thanks to the sym- 
metry of the problem, the ground state of Eq. (^) has 
a well defined positive parity. This allows us to elimi- 
nate redundant modes |^7| , saving memory and reaching 
higher energies and nonlinearities which otherwise would 
be computationally hard to attain. On the other hand, 
we have always checked that this method produced the 
same results as the complete one for a selected and sig- 
nificant set of parameter values. 

And finally it is important to note that the four-index 
tensor ( |l4| ) is indeed a product of two smaller tensors, 
corresponding to the integration on the z and r variables 
Jl6| . This decomposition is most important when work- 
ing with a large number of modes, because then the size 
of A becomes extremely large (i.e., 1600 4 elements for 
1600 modes). 

Concerning the evaluation of very high order polyno- 
mials as the ones involved in our computations it is nec- 
essary to say that it is not a simple task, specially for in- 
termediate values of the spatial variables since then there 
is a lot of comparable terms with usually different signs 
and the cancellations induce numerical instabilities. The 
usual procedure to avoid this difficulty is to use Horner's 
method p2|] to evaluate the polynomial, which is com- 
parable to using FFT techniques, but in our case this is 
not enough and the evaluation of higher order polynomi- 
als could be done only by the recursion formulas for the 
Hermite and Laguerre polynomials. 

We remark that the election of this spectral technique 
was largely influenced by the need of reaching high non- 
linearities which are not achievable using the other ap- 
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proaches. Further details on the numerical techni que as 
well as convergence proofs will be given elsewhere |24| 



C. Results for the stationary states and spectrum 

By using the preceding technique, we have searched 
the lowest states, (n z ,n r — 0), for each branch of the 
spectrum with a different vorticity, m = 0, . . . , 6. This 
was performed for two geometries corresponding to 7 = 
1 (spherically symmetric trap) and 7 = 2 (cigar shape 
trap), of a static trap, O = 0, while varying the intensity 
of the interaction from to approximately 50000. The 
results of this study are plotted in Fig. |l|. 



the lowest states, both the spectrum and the energies can 
be fitted to a simple formula 




FIG. 1. Plots (a) and (c) show the ground state energy, 
Eqo, and chemical potential, noo(U), dependence on the inter- 
action strength. Plots (b) and (d) show the chemical potential 
and the energy of the lowest state for each vorticity, always 
relative to value the ground state. The interaction values 
range from [7 = (upper diagonal) to U = 50000 (lowest di- 
agonal). All calculations shown correspond to the spherically 
symmetric trap, 7 = 1. 

Remarkably, in the absence of rotation, and up from 



Vom(N) = hoq(N) + u eS (N)m. 



(15) 



The first term is the chemical potential of the m = 
ground state and is not relevant to the dynamics. Using 
the Thomas-Fermi approximation one can show that it 
grows proportionally to /i cx N 2 ^ 5 , a behavior which is 
approximately reflected in the numerical results shown in 
Fig. 0(c). 

The second term is much more relevant to the evolu- 
tion of the condensate. It grows linearly, as the energy 
levels of a linear harmonic oscillator, with an effective 
frequency, w e //(7V), that decreases with the interaction. 
The fact that the higher levels of the spectrum of fi re- 
main equispaced even for large interactions is the reason 
why the condensate exhibits an exponentially divergent 
response to the parametric perturbation of the trap fre- 
quencies, as it is shown in Ref. |25| and |26| ]. 

Now we want to study the stationary solutions in the 
presence of rotation. For f2 ^ the proper functions with 
definite vorticity remain the same, while the chemical 
potential and the energy suffer a shift that depends on 
the vorticity of the state 



E nm (U,Q) = E nm (U,0)-ma 



(16) 




FIG. 2. Dependence of the energy levels on f2, Eo m (U, ft), 
for a fixed value of the interaction strength, U ~ 8000, and 
a spherically symmetric trap, 7 = 1. The horizontal line 
represents the vortex-free state, m — 0, the dashed line the 
m — 1 vortex state, and the dotted lines other multicharged 
vortex states. 

This shift gives rise to an ample phenomenology which 
is pictured in Fig. ^.First, we see that the degeneracy 
with respect to m is broken. The only other possible 
degeneracy that remains is with respect to the r and z 



4 



variables, but this will be removed for the case without 
spherical symmetry, 7^1. 

Second, the m = 1,2,3... branches of the spectrum 
become a minimum of the energy functional with respect 
to other branches for continuous intervals of the angular 
velocity, [Q m ,ft m+ i\, where 



n, 



(17) 



However this does not mean that in these intervals the 
771-th vortex state becomes a global minimum. Indeed, 
in Sect. Ill we will only be able to prove that only the 
77i = 1,2 vortex lines achieve the status of local minimum. 
It still remains an open question under which situations 
the ground state must have a well defined vorticity. 

Third, even though the separation between the m = 
and m = 1 states becomes very narrow for large interac- 
tions, the stabilization frequency Q.\ only approaches zero 
asymptotically with U. As a consequence, m = 1 states 
are never a global minimum of the energy in a stationary 
trap, a fact that can be checked by just inspecting the 
energy functional. 




FIG. 3. Frequency of stabilization for the vortex states 
in a spherically symmetric trap, 7 = 1, as a function of the 
nonlincarity. The lines are arranged in order of increasing 
vorticity, from m = 1 (solid line) to m = 6 (dashed line on 
the top). 

And finally, there is a critical value of for which 
the energy functional becomes unbounded by below [See 
Fig. |^] and which coincides with the separation between 
energy levels for large values of the vorticity. This critical 
value of the frequency, tt c , is such that all of the ground 
states for each value of the vorticity have the same energy 



Eom{U,Sl c ) =E ok (U,n c ), Vfc,m. 



(18) 



one finds that it is always smaller than the critical fre- 
quency of the linear case 



(19) 



III. STABILITY OF STATIONARY STATES 
A. The linear stability equations 

In the preceding section we obtained stationary so- 
lutions of the mean field model for the Bose-Einstein 
condensate, all of which had a well defined value of the 
third component of the angular momentum operator. We 
named those states vortices. In this section we now want 
to study the stability of these solutions according to sev- 
eral criteria of a local nature: local energetic stability 
and linear stability. 

We begin our study from the adimensionalized Gross- 
Pitaevskii equation ([}]). First we expand the condensate 
wavefunction around a stationary solution with a fixed 
vorticity. 

* (r, z, 0, t) =tp + eipi 

= [f(r, z)e imB + ea(r, z, 9, t)] e - *^*. (20) 

We insert this expansion in Eq. (||) and truncate the 
equations up to 0(e 1 ) thus getting 

id t a = [H + ittde + 2Uf 2 ] a + Uf 2 e - 2m9 a, (21a) 
-id t a = [Ho - ittde + 2Uf 2 } a + Uf 2 e 2m9 a, (21b) 

with Ho = -f A + i(7 2 r 2 + z 2 ) - fi(fi). We can also 
write this equation in a more compact form 



d - 

i—W = a z H(n)W = B(Q)W. 



by using the definitions 



W 



1 

-1 



W(fi) - H Q 



ind e + 2Uf 2 
Upe 2mB 



-iQd e + 2Uf 2 



(22) 

(23a) 
(23b) 
(23c) 



In the rest of this work we wish to study the dynamics 
that is involved in Eq. (|2^) . The simplest way to achieve 
this is to find a suitable basis in which Eq. (E3) becomes 
diagonal or almost diagonal. In other words, we want a 
set of vectors, W£ — (uk(r),Vk(r j) such that 



\ k W k = BW k . 



(24) 



Using Eqs. (|18|) and a fit similar to the one in Eq. (15) 



If B has such a diagonal Jordan form, then the pertur- 
bation evolves simply as 
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a(r,t) = ^c fc e iAt u fc (r,i). 



(25) 
(26) 



On the other hand, the lack of a diagonal form, or the 
existence of complex eigenvalues leads to instability in a 
way that we will precise later. 

Associated to Eq. (|2^) there is an energy functional, 

E 2 {a) = [ 2aH Q a + ^a 2 +^ 2 a 2 +4\^ \ 2 aa, (27) 



(28) 



and a constrained energy functional 

C 2 (a) = E 2 (a) - fi J \a\ 



which are the 0(e 2 ) terms in the expansion of (^) and 
(0), i.e. the energy introduced in the system by the per- 
turbation. If a diagonal Jordan form like the one of ( pij ) 
is possible, then it is easy to check that the second func- 
tional becomes diagonal, too 



£ 2 (a) =^M 2 A fe G(VK fe ), 
G(W k ) = [\u k \ 2 -\v k \ 2 . 



(29a) 
(29b) 



If the stationary state, i/jq, is a local minimum of the en- 
ergy subject to the constrain of a fixed norm (|^), then C 2 
must be positive for all perturbations, which has serious 
implications for the eigenvalues and eigenstates. We will 
refer to this later. 

When studying the condensate using tools from Quan- 
tum Field Theory, one may try a similar procedure [fl6|| , 
which is known as Bogoliubov's theory. In that frame- 
work, the a and a are linear operators in a Fock space, 
and one searches an expansion of these operators in terms 
of others that diagonalize the energy functional (|2^) and 
the evolution equations (|2^ ) . The resulting equations for 
the coefficients are known as Bog oliubov's equations and 
correspond to the equations (|24|) for u k and v k - 



B. Operational procedure 

It is now useful to perform an expansion of a and a into 
states of fixed vorticity so that the modes are separated 
into subspaces according to their vorticities 



TV, 



(30) 



These subspaces are not mixed by the action of the oper- 
ators of (^3|), and we can define their restriction to these 
subspaces 



B n (Q,)=a z n n (Q,), 
H n {O)=H 0n {O)+U n 



(31a) 
(31b) 



_(H»-(n-m)n \ , , 

n - ^ H 2m-n — {jyi — Tz)i7 ) ' ^ iCj 

H n = -\a + i( 7 V + z 2 ) + ^ + f MO) (31d) 



U n = Uf 



1 1 



(31e) 



With these definitions the diagonalization procedure 
Ml) becomes 



(n) _ 



= B n (Q)W, 



(n) 



n > m. 



(32) 



If G{W^ ] ) > 0, then 



(n) (2m— n) 



) is a Bogoliubov 

mode with energy e = A^™' and vorticities (n, 2m — 
n), whereas if G{W^) < then the excitation is 

(v^ m n \u^) with energy e = — Al . As a rule of 
thumb, the u function must always be the one with the 
largest contribution, which is formally stated in G(W) > 
0. In the following we will refer to these branches of the 
spectrum by the pairs of quantum numbers (n, 2m — n) 
and (2m — n, n), respectively. 

One may find, in principle, two kinds of solutions. 
First, the Bogoliubov operator may have complex eigen- 
values or even have a non diagonal Jordan form. In both 
cases we speak of dynamical instability because an arbi- 
trarily small perturbation departures from the original 
state exponentially or polynomially in time. Second, the 
linearized operator may have only real eigenvalues which 
should be interpreted as the change of energy in the con- 
densate due to excitations [See Eq. (psj)]. If A > the 
state under study t/'o is a local minimum of the energy 
functional (^) with respect to this family of perturba- 
tions, the A = case corresponds to the existence of 
degeneracy in the system, and finally if A < the sys- 
tem is told to be energetically unstable -i.e, excitations 
are energetically favorable and the state is not a local 
minimum of the energy. 

All of the five cases exposed above have the same impli- 
cations of stability for Eq. (^), which is a simple partial 
derivatives equation for an order parameter, and for the 
more complete Bogoliubov theory, where the perturba- 
tions are regarded as many-body corrections and involve 
more degrees of freedom. Nevertheless, it must be re- 
marked that of the two types of instability that can be 
found, i.e. dynamical and energetic instabilities, the sec- 
ond one is less harmful because it only affects the dynam- 
ics when there is some kind of dissipation that drives the 
system through the unstable branch. And even then the 
lifetime of the system can be significant if the intensity of 
the destabilizing mode is small compared to the typical 
times of evolution. 
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C. Numerical procedure 

We have discretized Eq. (|32|) in a basis which is es- 
sentially the same that we used to solve the stationary 
GPE. To be more precise, the expansion is as follows 



(n) 



k 



Pkr, 







Pl,2m-n 



(33) 



Here we have used the index convention explained above. 

In this basis, the operator H.Q n is diagonal, while the 
operator U can be calculated, either by means of integrals 
of the wavefunction itself in position representation, or 
by using a tensor of four indices which is similar to the 
one in Eq. ([hf). In any case, the equations are always 
linear, and so the study of the Bogoliubov spectrum con- 
sists in building and diagonalizing a large matrix of real 
numbers. 

Even though the procedure is quite simple, the ma- 
trices that one must build in order to resolve the case 
of strong interaction are very large and tend to exhaust 
computational resources. To be able to reach a large 
value of the nonlinearity we have had to work in a sub- 
space of states with even parity with respect to the Z 
axis. This way we could find the excitations with low- 
est energy for different vorticities, at the cost of missing 
those with odd parity, which are more energetic anyway 



D. Analytical results 

The study explained above does not have to be per- 
formed for all of the Bogoliubov operators in all of the 
possible situations. Here we will show several important 
results regarding when Eq. (^2|) may imply destabilizing 
modes. 

Lack of exponential instabilities in the Bogoliubov 
theory.- Any eigenvalue A satisfying Eq. (B2|) and 
G(W) ^ must be real. Eigenstates with G{W) = 
may involve complex eigenvalues but they are spurious 
and are introduced by the linearization procedure. 

This first part is shown simply by projecting the left 

and right hands of Eq. ([52]) against the vector w\ 
Omitting the indices the result is 



K")t 



Ar, 



\*-\vY) = \ (u n H n u + vH' m - n v) 



+ / Uf 2 \u + v\ 2 ~ (n 



mMM 2 + M 2 ) 



(34) 



The second part is more subtle. To prove it we must 
remember that solutions to Eq. (|2^) are stationary points 
of the action |p5[ , S = J L(t)dt corresponding to the 
following Lagrangian density 



L 



-{aott - aa t ) + C 2 {a). 



(35) 



Using (29a) it is easy to prove that the G(W) = modes 
are null modes that do not appear in the Lagrangian, and 
thus are not affected by the dynamics. 

It must be remarked that this result characterizes pos- 
sible eigenvalues, but does not grant that B n have a Bo- 
goliubov diagonalization. 

Sufficient condition for stability.- If the linearized 
Hamiltonian 7i n is positive definite, then B n may be di- 
agonalized, all of its eigenvalues are positive real numbers 
and there are no dynamical nor energetic instabilities. 

To prove this theorem one only needs to show that 
there's a one-to-one correspondence between the eigen- 

1/2 1/2 

functions of 7t n o~ z li. n and the eigenfunctions of o~ z TL n 
so that 



if and only if 



/2)| 



= \{K 



/C— 1/2) I 



.)) 



(36) 



(37) 



Then one uses this result to show that the eigenvalue in 
( |32| ) must be positive. 

Stability in stationary traps.- In Eq. (^2|), if ft = and 
n > 3m then the linearized Hamiltonian TL n is positive, 
the Bogoliubov operator B n can be diagonalized and it 
is also positive. Furthermore, if n > m then any real 
eigenvalue is positive, A > 0. 

The demonstration has four steps. First, one takes 
any value of n that satisfies that condition and proves 
that H 2m ' n > H m and H n > H m > 0. Second, this 
is used to prove that Ti.Q n > TLo m . Third, it shown that 
U n is positive which altogether implies TL n > 0. The last 
assertion may be easily checked with the help of (|34|) . 

The preceding two theorems imply that in a stationary 
trap any mode with negative energy must be comprised 
in the (0, —2m), . . . , (m, m) families, and any dynamic in- 
stability must lay in (0, —2m), . . . , (3m, 0). Thus we need 
only diagonalize a finite number of operators to make sure 
the system is stable or unstable. This result is an exten- 
sion of the one obtained in Ref. where a sufficient 
condition for stability is found to be n 2 > 4m 2 , without 
taking into account possible complex ei genv alues. 

Local stability under rotation. - In Eq. ( |32| ) the operator 
B n (fl) exhibits a linear dependence with respect to fl 



B n (ft) = B n (0) - (n-m)a 



(38) 



While the wave functions of the modes are the same as 
the ones of the stationary traps, the energies of the exci- 
tations suffer a global shift that depends on the vorticity 



A(ft) = A(0) - (n - m)0. 



(39) 



In general, the influence of these shifts has to be 
checked numerically. It is easy to show, however, that 
the shift is positive for n < m, which means that the pos- 
sibly negative eigenvalues in the range < n < m can be 
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suppressed if fl is large enough. Even more, as the shift 
is a real number, if one demonstrates that there are no 
dynamical instabilities in the stationary trap, then there 
will be no dynamical instabilities in the rotating trap, 
neither. 



E. Numerical results 

Summing up, from a practical point of view, the is- 
sue of stability consists in two different steps. The first 
one is the search for a stationary solution of the GPE 
with the appropriate vorticity, which we have already 
performed in Sect. |l[ and the second one is the study of 
the spectrum of the Bogoliubov operators for this partic- 
ular state. 

Stability of the m = l vortex-line in a stationary trap.- 
In this case of unit charge one only has to study a sin- 
gle operator, Bq, to know whether the system is stable. 
This calculation provides us with the branch of the spec- 
trum of excitations which is characterized by the quan- 
tum numbers (0,2) and (2,0), as already explained. We 
have done this for a wide range of nonlinearities in the ab- 
sence of rotation, £1 = 0, and the first conclusion is that 
the Bogoliubov operator has a diagonal Jordan form with 
all eigenvalues being real. 

In Fig U we show a selected set of the eigenvalues of 
the Bogoliubov operator, both for a spherically symmet- 
ric trap and an elongated trap. In those pictures one sees 
several things. First, there are two constant eigenvalues 
A = 1 which correspond to oscillations of the vortex line 
along the Z axis. Second, there is a single neutral mode 
A = only for the spherically symmetric trap, which cor- 
responds to the symmetry of rotation of the condensate 
around an axis on the XY plane. The symmetry and the 
mode disappear when 7 = 2 [See Fig. ||. And finally 
there is at least one negative eigenvalue fi < (more in 
the case of an elongated trap) which is responsible for 
the energetic destabilization of the system. The largest 
contribution to this destabilizing mode is a wavefunction 
captured in the vortex line and has zero vorticity (i.e. 
it is a core mode) [See Fig. as it was qualitatively 
predicted by Rokhsar in Ref. E| . 

We must remark that the number of unstable modes 
increases with the geometry factor: the more elongated 
the trap is, the easier is to transfer energy from the vor- 
tex to the core plus longitudinal excitations. In other 
words, for 7 < 1 (spherical or "pancake" traps) there's 
only one negative eigenvalue which corresponds to an ex- 
citation with a different vorticity than the unperturbed 
function, while for 7 > 1 we still have that mode, plus 
some more which are excited with respect to the Z axis. 
As a consequence, if the experiment is subject to dissi- 
pation and these unstable modes play a significant role 
in the dynamics, then the more elongated the trap is the 
less stable the vortex will be. 
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FIG. 4. Lowest eigenvalues of the Bogoliubov operator Bo 
for the m = 1 unperturbed state in (a) a spherically symmet- 
ric trap, 7=1, and (b) an axially symmetric trap, 7 = 2. The 
solid lines represent modes with quantum numbers (0, 2) and 
the dashed lines represent modes of the (2, 0) family. Crossing 
of levels is signaled with circles as a visual aid. 

In Fig ^| we also show the lowest eigenvalues of the 
families (-1, 3),(1, 1), (0,2), (2,0) and (-2,-4), that is, 
excitations where the main contribution is an eigenstate 
of L z with eigenvalues m — 0, ±1, ±2. In those pictures 
one sees that subspaces with excitations of the same vor- 
ticity but opposite sign have also different energy, a phe- 
nomenon which is solely due to the interaction. 
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which cannot be suppressed with any rotation below the 
critical value, f2 c > f2 > f^- 
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FIG. 5. Shape of the destabilizing mode. We show both 
the original solution (dotted line), the largest contribution, 
U, (dashed line) and the smallest contribution, v, (solid line). 
Functions have been rescaled to aid visualization. 

Stability of m = 1 vortex lines in rotating traps.- It 
was already proved ( |38| ) that, as the effect of rotation is 
gradually turned on, the modes with n < m and with 
ri > m are shifted up and down in the spectrum, respec- 
tively. It remained the question of whether the shift is 
enough to stabilize the vortex states, and the answer is 
yes, according to numerical experiments. 

First, as it is shown in Fig. |7|, the negative eigenvalue is 
slightly smaller than the stabilizing frequency, |Ao| < fii, 
which implies that for fl > J7i the energetically unstable 
branch with vorticity m = disappears. And second, 
the eigenvalues of B n for n > m are found to be larger 
than (n — m)VL\. In consequence for at least the inter- 
val [n^ffe) all of the B n operators are positive and the 
vortex with unit charge is a local minimum of the energy 
functional. 

In any case the shifts are always real, which implies 
that the B n operators remain diagonalizable with real 
eigenvalues and without dynamical instabilities. 

Stability of the m = 2 vortex line.- Another interesting 
configuration is the m = 2 multicharged vortex line. Here 
one suspects that a configuration with several vortices of 
unit charge has less energy than a single multicharged 
vortex, under all circumstances. In other words, they 
must be always energetically unstable. 

This intuitive perception is confirmed by the numerics. 
First the diagonalization of B\ reveals that this operator 
has at least one negative eigenvalue, while Bo has both 
negative eigenvalues and a pair of complex eigenvalues 
that, as we saw above, do not participate in the dynamics 
and must be discarded. Regarding the negative eigenval- 
ues, they do not decrease with the nonlinearity, but are 
always larger in absolute value than their linear limits. 
This implies that there are always negative eigenvalues 
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FIG. 6. Excitation energy of the lowest Bogoliubov modes 
for an unperturbed state with m — 1. The vorticities of each 
mode is written close to its corresponding line. Plot (a) cor- 
responds to 7 = 1, and (b) to 7 = 2. 

The immediate consequence of this linear stability 
analysis is that, due to the linearization of the energy 
(^7| ) not being positive, the m — 2 vortex-line is never a 
local minimum of the energy. This is true even for the 
parameter interval, [^2,^3), in which it has less energy 
than the rest of stationary states of well defined symme- 
try. If the m = 2 ground state is not a minimum, and 
the other symmetric states have more energy, we can con- 
clude that the minimum of the energy functional in the 
rotating trap with £ [f^, ^3) must be a state which is 
not symmetric with respect to rotations p9|. A similar 







analysis can be performed for the stationary states with 
m = 3, 4 . . . which extends this result to larger rotation 
frequencies, all below the critical one. 
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F. Lyapunov stability 

Speaking roughly, a solution of Eq. (||) is Lyapunov 
stable when every perturbed solution which is close 
enough to the original wave remains close throughout the 
evolution. The concepts of Lyapunov stability and lin- 
ear stability are close, but the latter does not imply the 
former as it is only defined in the limit of infinitesimal 
perturbations. 

Studying the Lyapunov stability of Eq. (||) theoret- 
ically is a difficult task that should be subject of fur- 
ther investigation. In the mean time we have performed 
an "empirical" study of the Lyapunov stability of the 
stationary solutions with m = 1 and m ~ 2 vortici- 
ties, by simulating numerically how they evolve for small 
perturbations and long times. The simulation was per- 
formed with a three-dimensional split-step pseudospec- 
tral method like the one from Ref . [^5| , using a 80 x 80 x 80 
points grid to study both the 7 = 1 and 7 = 2 problems. 

The main result of this complementary work is that 
both the unit charge vortex line and the multicharged 
vortex line are stable to perturbations which involve the 
destabilizing modes as defined by (|3^). For example, one 
may try to add a small contribution (0.5%) of a core 
mode to the m = 2 vortex, and with the result that 
the vortex line is split into two unit charge vortex lines, 
which rotate but remain close to the origin. We must re- 
mark that, although these simulations only work for finite 



times which are dictated by the precision of the scheme 
and the computational resources, these times are typi- 
cally 20 or 30 periods of the trap, which is much larger 
than any of the magnitudes that one may address theo- 
retically to the destabilization pro cess (i.e. the negative 
or complex eigenvalues of Eq. (|32|)) 

In the end, what this type of simulations reveal is that 
the m = 1 and m = 2 stationary states are energetically 
unstable, but this has no influence on the dynamic unless 
some other "mixing" or dissipative terms participate in 
the model. 



IV. CONCLUSIONS 

We have studied the vortex solutions of a dilute, 
nonuniform Bose condensed gas as modeled by the Gross- 
Pitaevskii equation (^) , both in a stationary, axially sym- 
metric trap, and subject to rotation (or a uniform mag- 
netic field). 

First, we have searched solutions of Eq. (||) that have 
the lowest energy and which are also eigenstates of the 
third component of the angular momentum operator, 
4>(r,z,8) = f(r,z)e" n0 , both in a stationary trap and 
in a rotating trap, and from small to very large nonlin- 
earities. It has been found that a nonzero angular speed 
(or magnetic field) is necessary in order to turn a vortex 
line state into a minimum of the energy functional with 
respect to other states of well-defined vorticity. However 
it remains open the question of whether the minimum of 
energy must have a well defined vorticity. 

Next we have studied the stability of these stationary 
solutions of the GPE. We have formulated a set of cou- 
pled equations that describe both the linearization of the 
GPE around a stationary solution and Bogoliubov's cor- 
rections to the mean field theory that describes the con- 
densate. It has been proved that the problem may not 
exhibit dynamical instabilities of exponential nature, plus 
several other theorems that describe the phenomenology 
associated to the possible instabilities. 

The perturbative equations have been solved numer- 
ically for stationary states having m = 1 and m = 2 
vorticities. In both cases it has been found that the only 
instability is of energetic nature, being limited to a small 
number of modes whose nature had already been pre- 
dicted in Q. 

For the vortex with unit charge we have found that 
this instability may be suppressed by rotating the trap 
at a suitable speed, and even when the trap is stationary, 
it is expected that it plays no significant role in the dy- 
namics unless there is enough dissipation as to take the 
system through the unstable branch. On the other hand, 
the linear stability analysis for the m — 2 multicharged 
vortex reveals that the energetic instability may never be 
suppressed, and that this configuration is never a min- 
imum of the energy functional, even though its lifetime 
is, once more, only conditioned by possible dissipation. 
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The last and probably most important conclusion of 
this work is that in the rotating trap, and for f2 > f^, 
the state of minimum energy is not an eigenstate of the 
L z operator, and thus it is not symmetric with respect 
to rotations. A similar result has been found in Ref. 
||l9f by means of a minimization procedure which is only 
justified in the limit of very small U, while our demon- 
stration remains valid for all nonlinearities, as far as the 
linearization procedure may be carried on. 

From a experimental point of view, this work has sev- 
eral implications. First, it is clear the conclusion that 
vortex lines with unit charge may be produced by ro- 
tating the trap at a suitable speed and then cooling the 
gas. Second, once rotation is removed, these vortices will 
survive for a long time if dissipation is small. Third, the 
multi-charged vortices are not minimum of the energy 
functional and thus it will be difficult to produce them 
by mean of cooling a rotating gas. And finally, if these 
multi-charged vortices are produced by some other mean 
such as Quantum Engineering, then we can assure that 
their lifetime will only depend on the intensity of dissi- 
pation, whose effect is to take the system either to the 
m = ground state if f2 < f2i, to the unit charge vortex- 
line state if Q < 0,2, or to a symmetry-less multicharged 
state if fi > 0,2 (A phenomenon which is regarded as 
splitting in the literature). 

The numerical results found in the paper have been 
possible due to the use of a powerful Galerkin spectral 
method optimized to allow the consideration of thou- 
sands of modes which is a step forward with respect to 
the previous analysis. 
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